The following script and analyses were written by Colleen Bove and James Umbanhowar. This version of the script was last committed to GitHub on 22 Febuary 2020

Results

In total, eleven studies met the standards of this meta-analysis, including the responses of thirteen Caribbean coral species collected from five different countries across the wider Caribbean (Figure 1 and Supplementary Table S1). Of the studies selected, only four performed fully factorial ocean acidification and warming experiments, and one performed two independent acidification and warming experiments. The most studied coral species from the Caribbean region were Porites astreoides (5 studies), Acropora cervicornis (4 studies), and Siderastrea siderea (4 studies). Finally, the Florida Keys and Belize were the two most-studied regions within the wider Caribbean fitting the criteria of this meta-analysis.


3.1 Overall calcification response

Meta-analysis of the dataset revealed that calcification rates of Caribbean corals were reduced by ocean warming but not ocean acidification (Figure 2 and Supplementary Figures S2). However, the 95% confidence interval of the combination of ocean warming and acidification overlapped zero, indicating no statistically clear trend towards synergistic or antagonistic effects of these treatments (Figure 2 and Supplementary Table S2).


Figure 1 Coral collection sites of all included studies with experimental study represented by colour and treatments represented by shape: acidification only (circle), warming only (triangle), and the combination of acidification and warming (square). The lower left insert displays close-up of the Belize collection sites and the upper right insert displays the Florida Keys collection sites.

Figure 1 Coral collection sites of all included studies with experimental study represented by colour and treatments represented by shape: acidification only (circle), warming only (triangle), and the combination of acidification and warming (square). The lower left insert displays close-up of the Belize collection sites and the upper right insert displays the Florida Keys collection sites.


# calculate ES and variance for ocean acidification
ES.oa <- escalc(measure = "SMDH", m1i = cm1, m2i = am3, sd1i = s1, sd2i = s3, 
    n1i = n1, n2i = n3, data = df)

# calculate ES and variance ES for ocean warming
ES.ow <- escalc(measure = "SMDH", m1i = cm1, m2i = wm4, sd1i = s1, sd2i = s4, 
    n1i = n1, n2i = n4, data = df)
ES.ow <- ES.ow %>% rename(yi.w = yi, vi.w = vi)

# create DF of calculated variables of OA
acid <- data.frame(cm1 = ES.oa$cm1, yi.a = ES.oa$yi, vi.a = ES.oa$vi)

# merge together warming and acid DF
ES.full <- merge(ES.ow, acid, by = "cm1")
ES.full <- ES.full[order(ES.full$cite), ]
ES.full$yi.w <- (ES.full$yi.w) * -1
ES.full$yi.a <- (ES.full$yi.a) * -1


### Standard deviation

# calculate sd for all
ES.full$sd1 <- ES.full$s1 * sqrt(ES.full$n1)
ES.full$sd3 <- ES.full$s3 * sqrt(ES.full$n3)
ES.full$sd4 <- ES.full$s4 * sqrt(ES.full$n4)
ES.full$sd6 <- ES.full$s6 * sqrt(ES.full$n6)

# calculate S for all treatments
ES.full$Sa <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd3^2/ES.full$n3))
ES.full$Sw <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd4^2/ES.full$n4))
ES.full$Si <- sqrt((ES.full$sd1^2/ES.full$n1) + (ES.full$sd6^2/ES.full$n6))


### Effect Size of Interaction term

# calculate ES using Harvey et al equations for acidification, warming, and
# interaction studies
ES.full$lnRRa <- ((ES.full$am3 - ES.full$cm1)/ES.full$Sa)
ES.full$lnRRt <- ((ES.full$wm4 - ES.full$cm1)/ES.full$Sw)
ES.full$yi.i <- ((ES.full$m6 - ES.full$wm4) - (ES.full$am3 - ES.full$cm1))/(2 * 
    ES.full$Si)

# calculate V for interaction studies ONLY
ES.full$vi.i <- (1/ES.full$n1) + (1/ES.full$n3) + (1/ES.full$n4) + (1/ES.full$n6) + 
    ((ES.full$yi.i^2)/(2 * (ES.full$n1 + ES.full$n3 + ES.full$n4 + ES.full$n6)))


write.csv(ES.full, file = "~/Bove2020_MetaAnalysis_FrontMarSci/Files/ES.full.csv")
# read in the saved CSV and then removed any unneeded columns
df<-read.csv("~/Bove2020_MetaAnalysis_FrontMarSci/Files/ES.full.csv")
df<-df[,-c(2, 12:13, 16:18, 20:22, 24:26, 41:49)] #remove any unneeded columns

### data manipulation

#subset df for variance (v) values only as new df
df2b<-subset(df, select=c(X,vi.i,vi.a,vi.w))

#subset df to remove only variance (v) values as another new df
df.m <-subset(df, select=-c(vi.i,vi.a,vi.w))

# this creates df with treat and V as column headers in long format
df2b<-gather(df2b, treat, v,vi.i:vi.w)

#substitute y for v for merge later
df2b$treat<-gsub('v','y', df2b$treat)

#creates another long format df with treat and Y as headers (of Y values)
df.m <-gather(df.m, treat, y,yi.w:yi.i)
df.m $treat<-as.factor(df.m $treat)

#merge two df into one
df2<-merge(df.m, df2b, by=c("X", "treat"))

#remove any NAs
#df2<-na.omit(df2)
df2$study <- paste(df2$author, df2$year, sep = "_")
df2$X<-NULL


## centering treatments

df2$mag.p <- df2$a.pco2 - df2$c.pco2
df2$mag.t <- df2$w.temp - df2$c.temp

df2$mag.cent.p <- df2$mag.p - 127 # make a centered around 0 column
df2$mag.cent.t <- df2$mag.t - 1.0 # make a centered around 0 column
# subset data by treatment
OW<-subset(df2, treat=="yi.w")
OA<-subset(df2, treat=="yi.a")
IN<-subset(df2, treat=="yi.i")

# fit model to each treatment
ow <- rma(y, v, data=OW)
oa <- rma(y, v, data=OA)
both <- rma(y, v, data=IN)

# wide to long format df
df2 <- gather(data = df2, key = "treat.mag", value = "mag", mag.p:mag.t)

# remove NA values from df
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

df2 <- completeFun(df2, "y")
df2 <- completeFun(df2, "mag")
global.mod <- rma.mv(y,v, mods = ~ treat , random = list(~ (treat|study), ~ (treat|spec)), data = df2) 

preds <- as.data.frame(predict(global.mod, newdata = df2))[,c(1,3:4)] # pull the mean and CI with the predict function
global.out <- cbind(df2, preds) # add columns to existing df for global output

#global.mod




global.out$treat <- factor(global.out$treat, levels = c("yi.a","yi.w","yi.i"))
global.out$wt <- 1/(global.out$v)

global_CI <- 
  ggplot()+
  theme(panel.grid.major=element_blank(), legend.position="none", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  ylab("Mean effect size (SMD)")+
  xlab("")+
  #ggtitle("Figure 2. Global standardized mean difference by treatment")+
  scale_x_discrete(labels=c("yi.a" = "acidification", "yi.w" = "warming", "yi.i" = "combination"))+
  geom_hline(yintercept=0, linetype=2, colour="black")+
  geom_point(aes(x = treat, y = y, size = wt),  colour="grey", data = global.out, alpha = 0.5, position=position_jitter(width=0.1))+
  #geom_errorbar(aes(ymin=lb, ymax=ub, width=0.1), size=0.6)+
  geom_linerange(data = global.out, aes(x = treat, ymin=ci.lb, ymax=ci.ub), size=0.8)+
  #ylim(-1.6,0.6)+
  geom_point(data = global.out, aes(x = treat, y = pred), size=10, shape=95)

div(ggplotly(global_CI), align = "center")
ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure2_GlobalCI.pdf", width = 5, height = 3.5)

Figure 2 Mean effect (standard mean difference) and 95% confidence interval of ocean acidification, warming, and the combination of acidification and warming on calcification rate for all studies in the meta-analysis. Grey circles indicate the effect size of each individual study and the size of each circle represents the weight of each study (1/SE). Clear statistical evidence of a treatment effect is identified when the 95% confidence interval does not overlap zero.



3.2 Calcification response of Florida Keys versus Belize corals

Corals from Belize only exhibited clearly reduced calcification rates under ocean warming (Figure 3A and Supplementary Table 3), while acidification, warming, and the combination of both stressors did not clearly alter experimental calcification rates of corals from the Florida Keys (Figure 3B and Supplementary Table S3). Further, the resulting QE suggests there is significant between-study variation (Supplementary Tables S3).


df2a<-subset(df2, region!="Curacao")
df2a<-subset(df2a, region!="Little Cayman")
df2a<-subset(df2a, region!="Mexico")

region.mod <- rma.mv(y,v, mods = ~ treat * region, random = list(~ (treat|study), ~ (treat|spec)), data = df2a)

preds2 <- as.data.frame(predict(region.mod, newdata = df2a))[,c(1,3:4)] # pull the mean and CI with the predict function
region.out <- cbind(df2a, preds2) # add columns to existing df for global output

#region.mod
region.out$treat <- factor(region.out$treat, levels = c("yi.a","yi.w","yi.i"))
region.out$wt <- 1/(region.out$v)

region_CI <- ggplot()+
  theme(panel.grid.major=element_blank(), legend.position="none", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  ylab("Mean effect size (SMD)")+
  xlab("")+
  scale_x_discrete(labels=c("yi.a" = "acidification", "yi.w" = "warming", "yi.i" = "combination"))+
  geom_hline(yintercept=0, linetype=2, colour="black")+
  geom_point(aes(x = treat, y = y, size = wt),  colour="grey", data = region.out, alpha = 0.5, position=position_jitter(width=0.1))+
  #ggtitle("Figure 3.Standardized mean difference by treatment for Florida Keys and Southern Belize")+
  #geom_errorbar(aes(ymin=lb, ymax=ub, width=0.1), size=0.6)+
  geom_linerange(data = region.out, aes(x = treat, ymin=ci.lb, ymax=ci.ub), size=0.8)+
  geom_point(data = region.out, aes(x = treat, y = pred), size=10, shape=95)+
  facet_wrap(~region)  

div(ggplotly(region_CI), align = "center")
ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure3_RegionalCI.pdf", width = 8, height = 3.5)

Figure 3 Mean effect (standard mean difference) and 95% confidence interval of ocean acidification, warming, and the combination of acidification and warming on calcification rate for (A) Florida Keys corals and (B) Belize corals. Grey circles indicate the effect size of each individual study and the size of each circle represents the weight of each study (1/SE). Clear statistical evidence of a treatment effect is identified when the 95% confidence interval does not overlap zero.



3.3 Temperature and aragnoite saturation state impacts on calcification rates across studies

Secondary analysis of mean calcification rates (mg cm−2 day−1) against treatment temperature across all Florida and Belize studies revealed a parabolic response to temperature (Figure 4 and Supplementary Tables S4, S5). Similarly, mean calcification rates across ΩArag resulted in a nonlinear response to acidification (Figure 5 and Supplementary Tables S6, S7). Both nonlinear trends in response to temperature and ΩArag were a result of treatment rather than region (Supplementary Tables S6, S7), suggesting regional differences identified in the meta-analysis were due to experimental designs employed to represent current regional environmental differences.

calc.w<-calc[,-c(15:18, 23:27)] #remove columns not needed for warming analysis

# remove and NA columns from warming column
calc.w <- completeFun(calc.w, "wm4")

# remove mexico studies
calc.w<-subset(calc.w, region != "Mexico")

# subset calc for calcication rates
calc.b<-subset(calc.w, select=c(num, cm1, wm4))
# create df with treat and rate as column headers in long format
calc.b<-gather(calc.b, treat, rate, cm1:wm4)
# rename for useful treatment names
calc.b$treat<-gsub('cm1','control', calc.b$treat)
calc.b$treat<-gsub('wm4','warm', calc.b$treat)

#subset calc for SE
calc.c<-subset(calc.w, select=c(num, s1, s4))
# create df with treat and SE as column headers in long format
calc.c<-gather(calc.c, treat, SE, s1:s4)
# rename for useful treatment names
calc.c$treat<-gsub('s1','control', calc.c$treat)
calc.c$treat<-gsub('s4','warm', calc.c$treat)

#subset calc to remove only columns from above
calc.d <-subset(calc.w, select=-c(s1 ,s4 ,cm1, wm4, c.pco2, n1, n4))
# this creates df with treat and V as column headers in long format
calc.d<-gather(calc.d, treat, temp, c.temp:w.temp)
# rename for useful treatment names
calc.d$treat<-gsub('c.temp','control', calc.d$treat)
calc.d$treat<-gsub('w.temp','warm', calc.d$treat)

#merge dfs into one, final warming
calc.m<-merge(calc.d, calc.b, by=c("num", "treat"))
calc.warm<-merge(calc.m, calc.c, by=c("num", "treat"))

# create column for weight
calc.warm$wt <- 1 / calc.warm$SE

# create study column
calc.warm$study <- paste(calc.warm$author, calc.warm$year, sep = "_")

# center temps
calc.warm$cent.temp <- calc.warm$temp - mean(calc.warm$temp)
calc.warm$cent.pco2 <- rep("NA", 60)
calc.warm$scale.temp <- scale(calc.warm$temp, center = TRUE, scale = TRUE)

warm.rate.region.mod <- lmer(rate ~ scale.temp + region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) #sing fit
warm.rate.region2.mod <- lmer(rate ~ scale.temp + region + I(scale.temp^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) 
warm.rate2.mod <- lmer(rate ~ scale.temp + I(scale.temp^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) 
warm.rate.mod <- lmer(rate ~ scale.temp + days+ (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) 
warm.region.mod <- lmer(rate ~ region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.warm,REML=F) #sing fit


#summary(warm.rate.region2.mod) # view summary output of model
#confint(warm.rate.region2.mod, method="boot")

warm.pred <- data.frame(temp=seq(26,32.01,length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100), pred=predict(warm.rate.region2.mod, newdata = data.frame(scale.temp=seq(min(calc.warm$scale.temp),max(calc.warm$scale.temp),length.out=100), days=seq(min(calc.warm$days),max(calc.warm$days),length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100)), re.form=NA)) 
global_OW_calc <-
ggplot(calc.warm, aes(x=temp, y=rate))+
  geom_point(aes(size=1/SE, shape=region, fill=region), colour="black")+
    theme(panel.grid.major=element_blank(), legend.position="right", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  scale_shape_manual("Region", values = c(21,22,24))+
  scale_fill_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  scale_color_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  geom_line(data = warm.pred, aes(y = pred, group = region, colour = region), size = 1)+
  ylab(bquote('Calcification rate (mg' ~cm^-2* ~day^-1*')'))+
  xlab(bquote('Temperature ('*degree*C*')'))
global_OW_calc

ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure4_CalcWarm.pdf", width = 7, height = 4)

Figure 4. Mean calcification rate (mg cm−2 day−1) of corals from each warming study by treatment temperature (°C) with the linear mixed effects model quadratic fit. Shape and colour of each point denotes study region (blue circle = Belize; brown square = Florida) and size of shape represents the weight (1/SE) of study.



calc.a<-calc[,-c(19:27)] #remove columns not needed for acid analysis
calc.a$p.mag <- calc.a$c.pco2 - calc.a$a.pco2 # pCO2 magnitude
calc.a <- completeFun(calc.a, "am3") # subset only for acidification studies 
calc.a <- subset(calc.a, region != "Curacao")


#subset calc for calcication rates
calc.b<-subset(calc.a, select=c(num, cm1, am3))
# create df with treat and rate as column headers in long format
calc.b<-gather(calc.b, treat, rate, cm1:am3)
# rename for useful treatment names
calc.b$treat<-gsub('cm1','control', calc.b$treat)
calc.b$treat<-gsub('am3','acid', calc.b$treat)

#subset calc for SE
calc.c<-subset(calc.a, select=c(num, s1, s3))
# create df with treat and SE as column headers in long format
calc.c<-gather(calc.c, treat, SE, s1:s3)
# rename for useful treatment names
calc.c$treat<-gsub('s1','control', calc.c$treat)
calc.c$treat<-gsub('s3','acid', calc.c$treat)

#subset calc to remove only columns from above
calc.d <-subset(calc.a, select=-c(s1 ,s3 ,cm1, am3, c.temp, n1, n3))
# this creates df with treat and V as column headers in long format
calc.d<-gather(calc.d, treat, acid, c.pco2:a.pco2)
# rename for useful treatment names
calc.d$treat<-gsub('c.pco2','control', calc.d$treat)
calc.d$treat<-gsub('a.pco2','acid', calc.d$treat)

#merge dfs into one, final acid
calc.m<-merge(calc.d, calc.b, by=c("num", "treat"))
calc.arag<-merge(calc.m, calc.c, by=c("num", "treat"))

# create column for weight
calc.arag$wt <- 1 / calc.arag$SE

# create study column
calc.arag$study <- paste(calc.arag$author, calc.arag$year, sep = "_")

# center temps
calc.arag$cent.pco2 <- calc.arag$acid - mean(calc.arag$acid)
#calc.arag$cent.pco2 <- rep("NA", 66)
calc.arag$scale.arag <- scale(calc.arag$acid, center = TRUE, scale = TRUE)

acid.rate.region.mod <- lmer(rate ~ scale.arag + region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) #sing fit
acid.rate.region2.mod <- lmer(rate ~ scale.arag + region + I(scale.arag^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) 
acid.rate2.mod <- lmer(rate ~ scale.arag + I(scale.arag^2) + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) 
acid.rate.mod <- lmer(rate ~ scale.arag + days+ (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) 
acid.region.mod <- lmer(rate ~ region + days + (1|study) + (1|spec), weights = 1/SE, data = calc.arag,REML=F) #sing fit

#summary(acid.rate.region2.mod) # continuing with this one for simplicity
#confint(acid.rate.region2.mod, method="boot")

acid.pred <- data.frame(acid=seq(0.7,5.8,length.out=100), region=rep(c("Belize", "Florida Keys"), length.out = 100),pred=predict(acid.rate.region2.mod, newdata = data.frame(scale.arag=seq(min(calc.arag$scale.arag),max(calc.arag$scale.arag),length.out=100), days=seq(min(calc.arag$days),max(calc.arag$days), length.out=100), region=rep(c("Belize", "Florida Keys"), length.out=100)),re.form=NA))
global_OA_calc <-
  ggplot(calc.arag, aes(x=acid, y=rate))+
  geom_point(aes(size=1/SE, shape=region, fill=region), colour="black")+
  theme(panel.grid.major=element_blank(), legend.position="right", panel.grid.minor=element_blank(), panel.background=element_blank())+
  theme(axis.line=element_line(color="black"))+
  theme(legend.key=element_blank())+
  scale_shape_manual("Region", values = c(21,22,24))+
  scale_fill_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  scale_color_manual("Region", values = c("#5ab4ac", "#d8b365", "#f5f5f5"))+
  geom_line(data = acid.pred, aes(y = pred, group = region, colour = region), size = 1)+
  ylab(bquote('Calcification rate (mg' ~cm^-2* ~day^-1*')'))+
  xlab(bquote('Aragonite Saturation State'))
global_OA_calc

ggsave("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Manuscript/Figure5_CalcAcid.pdf", width = 7, height = 4)

Figure 5. Mean calcification rate (mg cm−2 day−1) of corals from each acidification study by treatment aragonite saturation state (ΩArag) with the linear mixed effects model quadratic fit. Shape and colour of each point denotes study region (blue circle = Belize; brown square = Florida Keys) and size of shape represents the weight (1/SE) of study.



3.4 Experimental design impacts on coral calcification rate in studies

Quantification of experimental design parameters within warming studies identified that magnitude of treatment, irradiance, seawater type used, and feeding frequency all clearly impacted calcification rates (Table 1). Specifically, studies that utilized natural seawater and those with a larger difference between the control and treatment temperatures within a study exhibited higher effect sizes, suggesting a less negative effect of treatment. Studies that employed higher irradiance levels in their systems demonstrated more negative effects of treatment light level. Finally, studies that reported feeding their corals twice a week were less impacted by warming treatment than those feeding three times a week, however, studies with no data on feeding were the least affected by treatment. Duration of experiment was deemed redundant in the model and was thus dropped.

Within the acidification studies, irradiance, seawater type used, feeding frequency, duration, and the interaction of duration with treatment magnitude impacted effect sizes, while magnitude alone was not clearly different (Table 1). Studies using natural seawater, employing higher irradiance levels, and those with longer durations resulted in great effect sizes, suggesting they lessened the effects of acidification treatment on calcification responses. Similar to warming studies, acidification studies in which feeding corals was conducted twice a week exhibited less negative responses to treatment than those feeding three times a week, with studies reporting no feeding data exhibiting the least negative responses to treatment. Finally, coral calcification responses were less impacted by acidification in studies with longer duration of exposure and a greater pCO2 change.



Table 1. Effect size estimate and 95% confidence intervals of experimental design parameters calculated for all wamring and acidification experiments.

full.param <- rbind(w.global.out[,-4], a.global.out[,-4])

param.tab <- kable(full.param, digits = 3) %>% 
  kable_styling(font_size = 14, full_width = F) %>% 
  pack_rows("Warming Experiments", 1, 6) %>%
  pack_rows("Acidification Experiments", 7, 13)
param.tab
Estimate Lower 95% CI Upper 95% CI col
Warming Experiments
Intercept -29.392 -39.564 -19.221 negative
Magnitude 10.341 5.580 15.101 positive
Irradiance -0.237 -0.336 -0.137 negative
Seawater 29.767 19.719 39.815 positive
Feeding (3x) -3.544 -4.967 -2.120 negative
Feeding (N.D.) 87.160 49.685 124.635 positive
Acidification Experiments
Intercept -2.033 -2.785 -1.281 negative
Magnitude 0.000 -0.002 0.003 positive
Irradiance 0.009 0.004 0.014 positive
Seawater 2.392 1.144 3.640 positive
Feeding (3x) -15.042 -22.607 -7.477 negative
Feeding (N.D.) 3.611 0.280 6.942 positive
Duration 0.206 0.090 0.323 positive
Magnitude x Duration -0.001 -0.001 0.000 negative

Supplemental Materials


Supplemental Tables


Supplementary Table S1. List of included studies in meta-analysis with experimental design parameters. Treatment temperature and pCO2 reported are for treatment levels used in this analysis. Missing data from studies is denoted by N.D.

table1 <- kable(meta_tab2, booktabs = T) %>%
  kable_styling(font_size = 10, full_width = F) %>% 
  column_spec(5, italic = T)
table1
Study Location Collection.Sites Species Treatments Experiment.Duration Feeding Seawater Control.Temperature Temperature Control.pCO2 pCO2 irradiance
1 Albright et al. 2011 Florida Little Grecian, Key Largo Porites astreoides pCO2 49 d NA Natural seawater 28.0 NA 380 800.0 NA
2 Bedwell-Ivers et al. 2017 Little Cayman Little Cayman Research Center Acropora cervicornis Porites divaricata pCO2 28 d NA Natural seawater 30.0 NA 487 791.0 800
3 Bove et al. 2019 Belize Inshore, Offshore (Sapadilla Cayes) Porites astreoides, Pseudodiploria strigosa, Siderastrea siderea, Undaria tenuifolia Warming, pCO2, combination 93 d Frozen and newly hatched Artemia every other day Filtered non-native seawater 28.0 31.0 405 701.0 300
4 Castillo et al. 2014 Belize Nearshore, Backreef, Forereef (Sapadilla Cayes) Siderastrea siderea Warming, pCO2 95 d Frozen Artemia every other day Artificial seawater 28.0 32.0 472 704.0 250
5 Enochs et al. 2014 Florida Dade County Acropora cervicornis pCO2 42 d NA Natural seawater 28.0 NA 513 870.9 178
6 Grottoli et al. 2014 Mexico Puerto Morelos Porites divaricata, Porites astreoides, Orbicella faveolata Warming 15 d Fed Natural seawater 30.6 31.6 NA NA 600
7 Horvath et al. 2016 Belize Nearshore, Backreef, Forereef (Sapadilla Cayes) Siderastrea siderea Warming, pCO2, combination 60 d Frozen Artemia twice weekly Artificial seawater 28.0 31.0 426 890.0 250
8 Jury et al. 2010 Curacao Seaquarium Madracis auretenra pCO2 2 h Artemia nauplii twice weekly Natural seawater 28.0 NA 391 1480.0 200
9 Okazaki et al. 2017 Florida Key Biscayne, Biscayne National Park, Key West Acropora cervicornis, Agaricia agaricites, Dichocoenia stokesii, Montastraea cavernosa, Orbicella faveolata, Porites astreoides, Porites divaricata, Pseudodiploria strigosa, Siderastrea siderea, Siderastrea radians, Solenastrea hyades Warming, pCO2, combination 6 w Live rotifers and larval twice weekly Natural seawater 27.0 30.3 400 1340.0 327
11 Towle et al. 2015 Florida South Florida Acropora cervicornis Warming, pCO2, combination 8 w Dried zooplankton twice weekly Natural seawater 26.0 30.0 390 800.0 353
12 Kenkel et al. 2013 Florida Inshore and Offshore Sugarloaf Key Porites astreoides Warming 43 d NA Natural seawater 27.2 30.9 NA NA NA
#save_kable(table1, "Table1_Studies_16Jan20.pdf")


Supplementary Table S2. Global meta-analysis mixed effects model (function rma.mv) output by treatment only, with random effects of study and species, used in Figure 2. The test for residual heterogeneity and significance is represented by QE (P-value) and the test of moderators is QM (P-value).

tidy.rma <- function(x, ...) {
  with(
    x,
    data.frame(
      estimate = b,
      std.error = se,
      z.value = zval,
      p.value = pval,
      conf.low = ci.lb,
      conf.high = ci.ub
    )
  )
}
tidy.rma.Q <- function(x, ...) {
  with(
    x,
    data.frame(
      QE = QE,
      p.value = QEp,
      QM = QM,
       p.value = QMp
    )
  )
}

global.Q <- tidy.rma.Q(global.mod)
global.Q <- data.frame(t(global.Q))
row.names(global.Q) <- c("QE", "p value", "QM", " p value")
colnames(global.Q) <- "estimate"
global.Q$std.error <- rep("", 4)
global.Q$z.value <- rep("", 4)
global.Q$p.value <- rep("", 4)
global.Q$conf.low <- rep("", 4)
global.Q$conf.high <- rep("", 4)


global.mod.tab <- tidy.rma(global.mod)
#row.names(global.mod.tab) <- c("acidification", "combination", "warming")
global.mod.tab$std.error <- round(global.mod.tab$std.error, 3)
global.mod.tab$z.value <- round(global.mod.tab$z.value, 3)
global.mod.tab$p.value <- round(global.mod.tab$p.value, 3)
global.mod.tab$conf.low <- round(global.mod.tab$conf.low, 3)
global.mod.tab$conf.high <- round(global.mod.tab$conf.high, 3)
global.mod.tab <- rbind(global.mod.tab, global.Q)

global.mod.tab <- global.mod.tab[,1:4]

table2 <- kable(global.mod.tab, digits = 2) %>%
  kable_styling(font_size = 14, full_width = F) %>% 
  pack_rows("Model Results", 1, 3) %>%
  pack_rows("Variance Components", 4, 7)
table2
estimate std.error z.value p.value
Model Results
intrcpt -1.20 0.776 -1.545 0.122
treatyi.i 1.19 1.41 0.843 0.399
treatyi.w -1.07 1.19 -0.898 0.369
Variance Components
QE 823.95
p value 0.00
QM 2.40
p value 0.30
#save_kable(table2, "Table2_GlobalMod_16Jan20.pdf")


Supplementary Table S3. Meta-analysis mixed effects model (function rma.mv) output of treatment by region (Belize versus Florida Keys), with random effects of study and species, used in Figure 3. The test for residual heterogeneity and significance is represented by QE (P-value) and the test of moderators is QM (P-value).

regional.Q <- tidy.rma.Q(region.mod)
regional.Q <- data.frame(t(regional.Q))
row.names(regional.Q) <- c("QE", "p value", "QM", " p value")
colnames(regional.Q) <- "estimate"
regional.Q$std.error <- rep("", 4)
regional.Q$z.value <- rep("", 4)
regional.Q$p.value <- rep("", 4)
regional.Q$conf.low <- rep("", 4)
regional.Q$conf.high <- rep("", 4)


region.mod.tab <- tidy.rma(region.mod)
#row.names(region.mod.tab) <- c("acidification", "combination", "warming","acidification ", "combination ", "warming ")
region.mod.tab$std.error <- round(region.mod.tab$std.error, 3)
region.mod.tab$z.value <- round(region.mod.tab$z.value, 3)
region.mod.tab$p.value <- round(region.mod.tab$p.value, 3)
region.mod.tab$conf.low <- round(region.mod.tab$conf.low, 3)
region.mod.tab$conf.high <- round(region.mod.tab$conf.high, 3)
region.mod.tab <- rbind(region.mod.tab, regional.Q)

region.mod.tab <- region.mod.tab[,1:4]

table3 <- kable(region.mod.tab, digits = 4, booktabs = T) %>%
  kable_styling(font_size = 14, full_width = F) %>% 
  pack_rows("Belize", 1, 3) %>%
  pack_rows("Florida", 4, 6) %>%
  pack_rows("Variance Components", 7, 10)
table3
estimate std.error z.value p.value
Belize
intrcpt -0.9288 1.324 -0.701 0.483
treatyi.i 1.5625 2.24 0.697 0.486
treatyi.w -3.6439 2.069 -1.761 0.078
Florida
regionFlorida Keys 0.0861 1.879 0.046 0.963
treatyi.i:regionFlorida Keys -0.9606 3.184 -0.302 0.763
treatyi.w:regionFlorida Keys 5.1661 3.061 1.688 0.091
Variance Components
QE 504.3862
p value 0.0000
QM 9.3736
p value 0.0951
#save_kable(table3, "Table3_RegionalMod_16Jan20.pdf")


Supplementary Table S4. Temperature and region linear mixed effects model section using AICc. All models were run with random effects for study and species and a weight of 1/SE. The asterisk (*) denotes the selected model run for the final analysis.

warm.mag.tab <- AICc(warm.rate.region.mod, warm.rate.region2.mod, warm.rate2.mod, warm.rate.mod, warm.region.mod)

tableS4 <- kable(warm.mag.tab, digits = 2, booktabs = T) %>%
  kable_styling(font_size = 14, full_width = F)  %>% 
  row_spec(2, bold = T)
tableS4
df AICc
warm.rate.region.mod 7 109.99
warm.rate.region2.mod 8 75.39
warm.rate2.mod 7 72.85
warm.rate.mod 6 107.84
warm.region.mod 6 125.18


Supplementary Table S5. Temperature and region best fit linear mixed effects model output and 95% confidence intervals of the calcification rates in response to treatment temperature plotted in Figure 4.

tab_model(warm.rate.region2.mod, show.r2 = FALSE, show.aicc = TRUE)
  rate
Predictors Estimates CI p
(Intercept) 0.74 -1.72 – 3.21 0.555
scale.temp -0.23 -0.31 – -0.16 <0.001
region: Florida Keys 0.24 -1.05 – 1.53 0.713
I(scale.temp^2) -0.36 -0.44 – -0.27 <0.001
days 0.01 -0.02 – 0.04 0.659
Random Effects
σ2 0.28
τ00 spec 0.20
τ00 study 0.18
ICC 0.57
N study 5
N spec 11
Observations 60
AIC 75.386


Supplementary Table S6. Aragonite saturation state and region linear mixed effects model section using AICc. All models were run with random effects for study and species and a weight of 1/SE. The asterisk (*) denotes the bet fit model.

arag.mag.tab <- AICc(acid.rate.region.mod, acid.rate.region2.mod, acid.rate2.mod, acid.rate.mod, acid.region.mod)

tableS6 <- kable(arag.mag.tab, digits = 2, booktabs = T) %>%
  kable_styling(font_size = 14, full_width = F)  %>% 
  row_spec(2, bold = T)
tableS6
df AICc
acid.rate.region.mod 7 57.97
acid.rate.region2.mod 8 49.85
acid.rate2.mod 7 47.98
acid.rate.mod 6 56.39
acid.region.mod 6 105.28


Supplementary Table S7. Aragonite saturation state (ΩArag) and region best fit linear mixed effects model output and 95% confidence intervals of the calcification rates in response to treatment ΩArag plotted in Figure 5.

tab_model(acid.rate.region2.mod, show.r2 = FALSE, show.aicc = TRUE)
  rate
Predictors Estimates CI p
(Intercept) -0.19 -3.06 – 2.69 0.900
scale.arag 0.10 0.07 – 0.14 <0.001
region: Florida Keys 0.47 -1.00 – 1.95 0.532
I(scale.arag^2) -0.04 -0.07 – -0.02 0.001
days 0.02 -0.02 – 0.05 0.376
Random Effects
σ2 0.20
τ00 spec 0.17
τ00 study 0.25
ICC 0.68
N study 5
N spec 12
Observations 140
AIC 49.853



Supplemental Figures

oa.for <- forest(oa, slab=paste(OA$cite, OA$spec, OA$site, sep=" - "), order=order(OA$cite))

Supplementary Figure S1. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of ocean acidification only studies. Studies with confidence intervals not overlapping zero (dotted line) denote either a significant increase (positive values) or decrease (negative values) in calcification response under experimental acidification compared to the corresponding control treatment.

oa.fun <- funnel(oa, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)

Supplementary Figure S2. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of ocean acidification only studies. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.


pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS1_OAforest.pdf", width = 11, height = 12)
forest(oa, slab=paste(OA$cite, OA$spec, OA$site, sep=" - "), order=order(OA$cite))
dev.off()

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS2_OAfunnel.pdf", width = 7, height = 6)
funnel(oa, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
dev.off()
ow.for <- forest(ow, slab=paste(OW$cite, OW$spec, OW$site, sep=" - "), order=order(OW$cite))

Supplementary Figure S3. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of ocean warming only studies. Studies with confidence intervals not overlapping zero (dotted line) denote either a significant increase (positive values) or decrease (negative values) in calcification response under experimental warming compared to the corresponding control treatment.

ow.fun <- funnel(ow, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)

Supplementary Figure S4. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of ocean warming only studies. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.


pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS3_OWforest.pdf", width = 9, height = 10)
forest(ow, slab=paste(OW$cite, OW$spec, OW$site, sep=" - "), order=order(OW$cite))
dev.off()

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS4_OWfunnel.pdf", width = 7, height = 5)
funnel(ow, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend="topleft")
dev.off()
com.for <- forest(both, slab=paste(IN$cite, IN$spec, IN$site, sep=" - "), order=order(IN$cite))

Supplementary Figure S5. Forest plot depicting individual standard mean difference (SMD) and confidence interval of each species-site-study combination included in the meta-analysis of the combination of ocean acidification and warming studies.

com.fun <- funnel(both, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)

Supplementary Figure S6. Funnel plot of the standard mean difference (SMD) against standard error of each species-site-study combination included in the meta-analysis of combined ocean acidification and warming only. Dotted lines represent confidence intervals of all studies, background plot colour represents statistical significance of individual study, the black vertical line denotes the overall effect of all included ocean acidification studies.


pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS5_bothforest.pdf", width = 9, height = 8)
forest(both, slab=paste(IN$cite, IN$spec, IN$site, sep=" - "), order=order(IN$cite))
dev.off()

pdf("~/Bove2020_MetaAnalysis_FrontMarSci/Figures/Supplemental/FigureS6_bothfunnel.pdf", width = 7, height = 6)
funnel(both, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"), legend=TRUE)
dev.off()
sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.0.0      kableExtra_1.1.0   sjlabelled_1.1.1  
##  [4] sjmisc_2.8.2       sjPlot_2.7.2       plotly_4.9.0      
##  [7] shiny_1.4.0        Rmisc_1.5          lattice_0.20-38   
## [10] MuMIn_1.43.6       RColorBrewer_1.1-2 tidystats_0.4     
## [13] forcats_0.4.0      stringr_1.4.0      dplyr_0.8.3       
## [16] purrr_0.3.3        readr_1.3.1        tibble_2.1.3      
## [19] tidyverse_1.2.1    lme4_1.1-21        plyr_1.8.4        
## [22] tidyr_1.0.0        ggplot2_3.2.1      MAd_0.8-2.1       
## [25] metafor_2.1-0      Matrix_1.2-17     
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.0-10    minqa_1.2.4       colorspace_1.4-1 
##  [4] ellipsis_0.3.0    estimability_1.3  parameters_0.2.0 
##  [7] rstudioapi_0.10   ggrepel_0.8.1     mvtnorm_1.0-11   
## [10] lubridate_1.7.4   xml2_1.2.2        codetools_0.2-16 
## [13] splines_3.5.2     mnormt_1.5-5      knitr_1.25       
## [16] zeallot_0.1.0     jsonlite_1.6      nloptr_1.2.1     
## [19] ggeffects_0.12.0  Cairo_1.5-10      broom_0.5.2      
## [22] compiler_3.5.2    httr_1.4.1        sjstats_0.17.6   
## [25] emmeans_1.4.2     backports_1.1.5   assertthat_0.2.1 
## [28] fastmap_1.0.1     lazyeval_0.2.2    cli_1.1.0        
## [31] later_1.0.0       formatR_1.7       htmltools_0.4.0  
## [34] tools_3.5.2       coda_0.19-3       gtable_0.3.0     
## [37] glue_1.3.1        Rcpp_1.0.2        cellranger_1.1.0 
## [40] vctrs_0.2.0       nlme_3.1-141      crosstalk_1.0.0  
## [43] psych_1.8.12      insight_0.6.0     xfun_0.10        
## [46] rvest_0.3.4       mime_0.7          lifecycle_0.1.0  
## [49] MASS_7.3-51.4     zoo_1.8-6         scales_1.0.0     
## [52] hms_0.5.1         promises_1.1.0    parallel_3.5.2   
## [55] sandwich_2.5-1    yaml_2.2.0        stringi_1.4.3    
## [58] highr_0.8         bayestestR_0.4.0  boot_1.3-23      
## [61] rlang_0.4.1       pkgconfig_2.0.3   evaluate_0.14    
## [64] htmlwidgets_1.5.1 labeling_0.3      tidyselect_0.2.5 
## [67] magrittr_1.5      R6_2.4.0          generics_0.0.2   
## [70] multcomp_1.4-10   pillar_1.4.2      haven_2.1.1      
## [73] foreign_0.8-72    withr_2.1.2       survival_2.44-1.1
## [76] performance_0.4.0 modelr_0.1.5      crayon_1.3.4     
## [79] rmarkdown_1.16    grid_3.5.2        readxl_1.3.1     
## [82] data.table_1.12.6 digest_0.6.22     webshot_0.5.1    
## [85] xtable_1.8-4      httpuv_1.5.2      stats4_3.5.2     
## [88] munsell_0.5.0     viridisLite_0.3.0